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Abstract 

We show that the exact integrator for the classical Kepler mo- 
tion, recently found by Kozlov (J. Phys. A: Math. Theor. 40 (2007) 
4529-4539), can be derived in a simple natural way (using well known 
exact discretization of the harmonic oscillator). We also turn atten- 
tion on important earlier references, where the exact discretization of 
the 4-dimensional isotropic harmonic oscillator has been applied to 
the perturbed Kepler problem. 
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In recent years there were proposed several conservative discretizations 
of the classical Kepler problem [U EJ El lU E]- These numerical integrators 
preserve all integrals of motion and trajectories but only Kozlov's schemes 
are of order higher than 2. Kozlov found also the exact integrator by guessing 
its proper form and summing up some infinite series [I]. 

In this comment we show that Kozlov's exact integrator can be derived in 
a simple elementary way. Conservative discretizations of the 3-dimensional 
Kepler motion obtained in [H El H] consist in applying the midpoint rule 
(or the discrete gradient method, compare [6]) to the isotropic 4-dimensional 
harmonic oscillator equations: 

f = i p - f = 2£Q ' (Q,P e R 4 ) (i) 

where E = const is the energy integral of the considered Kepler motion. 
Then, the Authors of [31 0] use the Kustaanheimo-Stiefel (KS) transforma- 
tion. This classical transformation is given by (£7], see also [1]): 

Qi ~ Q\ ~ Q\ + Ql 

q = | 2Q1Q2 - 2Q 3 Q 4 I , (2) 
2Q 1 Q 3 + 2Q 2 Q 4 

P1Q1 ~ P2Q2 - P3Q3 + PaQa 

P1Q2 + P2Q1 - P3Q4 - P4Q3 I , (3) 
P1Q3 + P2Q4 + P3Q1 + P*Q2 
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2|Q| 

where Q, P are subject to the constraint 

P1Q4 - P2Q3 + P3Q2 - P4Q1 = . (4) 
The KS transformation, together with the Levi-Civita time transformation 

maps the 4-dimensional harmonic oscillator ([1]) into the 3-dimensional Kepler 
problem equations: 

rf q dp kq 3 

^ = P ' -dl = -W (q ' pGR) ' (6) 
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where k = const. Using (j2J), (jSJ) and (jlj) we can verify useful identities 

|q| 2 = |Q|\ |P| 2 = 4| P | 2 |Q| 2 , (7) 
which imply the equivalence of the energy conservation laws: 

\v 2 ~^ = E ^ i|P| 2 -£|Q| 2 = A;. (8) 

/ qj o 

The phenomenon of interchanching coupling constants with integrals of 
motion (like k E) is quite well known in the theory of integrable systems, 
see [8] (compare also [9], where more general results can be found). 

In order to derive Kozlov's numerical results in a simple straightforward 
way it is sufficient to notice that the KS transformation (used by Kozlov), 
reduces the Kepler motion to linear ordinary differential equations with con- 
stant coefficients (namely, to the harmonic oscillator) and for all such equa- 
tions there exist explicit exact numerical integrators ([I0JEET], see also [T2]). 
By the exact discretization of an ordinary differential equation x = f(x), 
where x(t) G M. N , we mean the one-step numerical scheme of the form 
X n+ i = $h(X n ), such that X n = x(t n ), compare [TUl fTT]. 

The system (JTJ, equivalent to the 4-dimensional harmonic oscillator equa- 
tion, admits the exact discretization (see, for instance, [12]): 

Qj+i - Qj = i Pj+i + Pj 

5(hj) 4 2 

(9) 

Pj+i ~ Pj 9F Qj+i + Qj 

5(M 2 

where /ij := Sj + \ — Sj is the (variable) s-step, Qj,Pj denote jth iteration of 
the numerical scheme (not to be confused with coordinates Qj, Pj), and 

^■) = -tan^, ^ = -\ E . (10) 

In the case of the constant step hj = h and E < 0, we recognize here the exact 
integrator found by Kozlov (see formulae (4.11) and (4.14) from [4j, taking 
into account that 5(h) = ha(h) = hb(h) and E = —A). The hyperbolic 
and parabolic cases (formulae (4.16) and (4.18) from [4]) follow immediately 
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when we take imaginary uj (i.e., E > 0) or u = 0, respectively. The exact 
numerical scheme preserves the energy integral, i.e., 



IP-I 2 



E\Qj\ 



k 



Note that the system ([9]) can be rewritten in the explicit form: 



Qj+i = cos lo hj Qj + 



sin uhj 



Au J ' 
Au sin to hj Qj + cos uihj Pj . 



This system is a direct consequence of evaluating the exact solution of 
s = Sj and s = Sj + hj, compare [i~2] . 



(12) 
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The equation (jSJ) can be solved exactly in different (but more or less 
equivalent) ways, compare [U [131 CH] . Here we propose one more approach, 
reducing this problem to linear ordinary differential equations with constant 
coefficients. If Q,P satisfy (CQ) and t satisfies (J5J), then we easily check that 



dw 

ds 



fhv, 



w 



/ IQI 2 \ 
|p| 2 

Q P 
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n 
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(13) 
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In such case we can proceed in a standard way. The general solution is given 
by w(s) = exp(sfi)w(0). Therefore, the exact discretization, w n = w(/m), 
satisfies 



w n+ i = exp(M7)w r 



(14) 



and the problem reduces to the well known, purely algebraic procedure of 
computing e nh . In our particular case we observe that f2 4 = 2EQ 2 which 
simplifies computations. The last row in the equation (I14p reads 



tj+i 



sin 2huj 
4u 



/"v 1 2 I j I | ^ ( I v 1 2 I jf I 

Qjl ~16^J + 2 U Qil + 16^ 



Q,- • P,- sin 2 hu 



Au 2 



(15) 
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One can check by direct computation that the discretization (1151) . although 
have a simpler form, is identical with the formulae (4.11), (4.15) of [4]. Fi- 
nally, eliminating |Pj| 2 by virtue of ffTTl) . we get 

hk ( sin2/iuA sm.2hu,^ l9 Q,- • P, sin 2 luo , , 

<*« = t ' + ^[ 1 -^r) + + ^-±7* — • < 16 > 

Another approach (see [13]) consists in computing the integral f \Q(s)\ 2 ds, 
where Q is the exact solution of ([JJ). The formula (86) from [13] is identical 
to ffTBI) (although notation is quite different). 

In celestial mechanics the exact discretization of the Kepler motion via the 
KS transformation appeared as a quite natural step [121 [111 [15] , although the 
conservative properties of the exact integrator were not discussed explicitly in 
these papers. A long time ago Stiefel and Bettis, working in the framwork of 
the Gautschi approach [17J, applied the exact discretization of the harmonic 
oscillator to the perturbed Kepler motion [T5J dS] . 

More recently, Mikkola [13J and Breiter [14] proposed new integrators for 
the perturbed Kepler problem, using the exact solution of the 4-dimensional 
isotropic harmonic oscillator equation and the exact discretization ffTBI) of 
the time (known to Stumpff even before the KS transform was invented, 
compare [13])- In particular, the numerical scheme (fT2"j) can be found in [13J, 
p. 162, and in [14], p. 234. Breiter follows [18] using an additional constant in 
the definition of the KS transformation (in fact scaling both q and p). The 
freedom of choosing this parameter can be used to fix the value of uj (e.g., 
uj = 1) which may have numerical advantages [T4] . 

These important results of celestial mechanics are not mentioned in [4] 
and, in general, they seem to be rather unknown in the field of geomet- 
ric numerical integration [19] . It is worthwhile to mention that the exact 
discretization of the harmonic oscillator equation has been recently used to 
construct new geometric integrators of high accuracy ( "locally exact discrete 
gradient schemes") [20] ■ We plan to apply such scheme to the perturbed 
Kepler problem using the Kustaanheimo-Stiefel map. 
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